Bi-stability in voltage-biased NISIN structures 



OO 
O 

o 

(N 

< 

(N 

o 
o 



c 

o 
o 



> 

OO 

in 
cn 

o 

OO 

o 



X 



I. Snyman 1 ' 2 and Yu. V. Nazarov 3 

1 National Institute for Theoretical Physics, Private Bag XI, 7602 Matieland, South Africa 
2 Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands 
s Kavli Institute of Nanoscience, Delft University of Technology, 2628 CJ Delft, The Netherlands 

(Dated: August 2008) 

As a generic example of a voltage-driven superconducting structure we study a short supercon- 
ductor connected to normal leads by means of low transparency tunnel junctions, with a voltage bias 
V between the leads. The superconducting order parameter A is to be determined self-consistently. 
We study the stationary states of the system as well as the dynamics after a perturbation. We 
find a region in parameter space where there are two stable stationary states at a given voltage. 
These bi-stable states are distinguished by distinct values of the superconducting order parameter 
A and of the current between the leads. We have evaluated (1) the multi-valued superconducting 
order parameter A at given V; (2) the current between the leads at a given V; and (3) the critical 
voltage at which superconductivity in the island ceases. With regards to dynamics, we find numer- 
ical evidence that the stationary states are stable and that no complicated non-stationary regime 
can be induced by changing the voltage. This result is somewhat unexpected and by no means 
trivial, given the fact that the system is driven out of equilibrium. The response to a change in the 
voltage is always gradual, even in the regime where changing the interaction strength induces rapid 
anharmonic oscillations of the order parameter. 
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I. INTRODUCTION 

Electron transport devices combining superconducting 
(S), insulating (I) and normal metal (N) elements are 
known as superconducting hetero-structures. Often such 
hetero-structures are more than the sum of their partsi^ 
Phenomena that are not present in bulk S, I or N systems 
appear when a device contains junction between these 
components. The following examples are well known: 
(1) The conductance of a high transparency NS junction 
does not equal the conductance of the normal metal on 
its own, as one might naively expect. If the normal metal 
is free of impurities, the conductance is higher than that 
of the normal metali 3 - This surprising effect is due to a 
process known as Andreev reflection^ During Andreev 
reflection at an NS interface, an electron impinging on 
the interface from the N side is reflected back as a hole, 
while a Cooper pair propagates away from the interface 
on the S side. (2) In Josephson junctions, the simplest of 
which is perhaps the SIS hetero-structure^ a DC current 
can flow at zero bias voltage. This happens when the 
superconducting phase difference across the junction is 
non-zero £ 

The above examples can be understood in terms of 
equilibrium properties of the hetero-structure. When a 
superconducting device is perturbed outside equilibrium, 
yet more interesting effects can occur) 7 - for instance, oscil- 
lations under stationary non-equilibrium conditions. An 
elementary example: if a Josephson junction is biased 
with a DC (i.e. fixed) voltage, an AC (i.e. oscillating) 
current flows through the junction^ Another example 
of the kind has been investigated in the context of cold 
Fermi gases in optical traps. In these systems, the in- 
teraction between atoms can be tuned and changed by 



means of a so-called Feshbach resonance. If the interac- 
tion is attractive, the gas forms a BCS-condensate. Re- 
cent studies^ have considered what happens if the value 
of the attractive pairing interaction is changed abruptly. 
It was discovered that, depending on the ratio between 
the initial and final values of the interaction strength, the 
condensate order parameter can perform anharmonic os- 
cillations that do not decay in time. 

The initial motivation for the research presented in this 
paper came form the study of Keizer et al.^ where the 
authors investigated the suppression of the superconduct- 
ing order parameter by a voltage applied to a supercon- 
ducting wire. It was assumed that A remains stationary. 
However, this assumption does not seem well-justified: 
the stationary voltage could induce periodic oscillations 
of |A| or even richer chaotic dynamics. Thus prompted, 
we wanted to address the validity of this assumption for a 
decidedly simpler NISIN structure, namely a short super- 
conductor connected to normal leads by means of tunnel 
junctions. The structure is biased with a voltage V. 

We require that (1) the dominant energy relaxation 
mechanism in the superconductor is the tunneling of 
electrons to the leads, and (2) spatial variations of the 
superconducting order parameter inside the supercon- 
ductor are negligible. To meet the first requirement, 
the superconductor must have dimensions smaller than 
the inelastic scattering length of quasi-particles. This 
is not an unrealistic requirement given current experi- 
mental techniques. To meet the second requirement, the 
superconductor should firstly contain impurities or have 
an irregular shape, so that the electron wave-functions 
of the isolated island are isotropic on the scale of the 
superconducting coherence length*^ Secondly, the tun- 
nel junctions connecting it to the leads should have a 
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bigger normal-state resistance than that of the supercon- 
ductor proper. In this case, opening up the system by 
connecting leads does not re- introduce spatial anisotropy 
of wave-functions inside the island. 

The study of NISIN structures has a long history 12 ' 13 
Our study complements several previous studies i 10 ' 14 ' 15 
These dealt with quasi-one-dimensional superconducting 
wires between normal leads. Setups where either the su- 
perconductor was impurity-free or the transparency of 
the NS interfaces were high were considered. For these 
setups, spatial variations of the order parameter, specif- 
ically the spatial gradient of the superconducting phase, 
can be large. Including these spatial variations in the 
description of the superconductor significantly compli- 
cates matters. Hence these studies focused on numeri- 
cal calculations and assumed that the superconducting 
order-parameter and all other quantities of interest were 
stationary. It should also be mentioned that asymmet- 
ric couplings, where the superconductor is coupled more 
strongly to one lead than the other, did not receive de- 
tailed analysis. The only asymmetric setup considered 
consisted of one interface with tunable transparency and 
the other perfectly transparent.— One of the main con- 
clusions of these studies is that, if the bias voltage is large 
enough, the system switches to the normal state. Some 
evidence for a bi-stable region where, depending on the 
history of the system, either the superconducting or the 
normal state can occur at a given voltage, was reported^ 

The absence of spatial variations in the system we 
study allows us to perform analytical calculations, pro- 
vided we assume stationarity. Results are obtained for an 
arbitrary ratio of the coupling strengths to the leads. We 
derive transcendental equations relating the supercon- 
ducting order parameter to the bias voltage, and derive 
an explicit formula for the current between the leads. As 
mentioned, the assumption of stationarity is however not 
a priori justified. As was seen in the examples mentioned 
at the beginning of this introduction, non-equilibrium 
conditions in superconductors often go hand in hand with 
non-stationary behavior of observable quantities. Indeed, 
the NISIN junction that we study is a non-linear system 
subjected to a driving force (and to damping). Non- 
linearity here means that the dynamical equations for 
one-particle Green functions are not linear in the Green 
functions. This is due to the existence of a non-zero 
superconducting order parameter. The driving force is 
provided by the voltage (and the damping by tunneling 
of electrons from the island into the leads). Non- linear 
driven systems (think of the nonlinear pendulum) often 
have chaotic dynamics. The assumption of stationarity 
would miss this. We therefore supplement our analyti- 
cal calculation with numerical calculations that study the 
dynamics in real time. 

Our main results are the following: The stationary 
states that we found analytically are stable. Further- 
more, there is a parameter region where two different 
stationary states are stable at the same voltage. (This is 
the "bi-stability" of the title.) For a symmetric coupling 



to the source and drain leads, one of the two states is 
superconducting (characterized by a non-zero order pa- 
rameter) and the other is normal. Since we are in the 
regime of high tunnel barriers, at a given voltage, the su- 
perconducting island allows less current to flow between 
the leads than the island in the normal stated This cur- 
rent is a directly measurable quantity and allows one to 
distinguish between superconducting and normal states. 
For some asymmetric couplings however, both the sta- 
ble states are superconducting. We have calculated the 
current that flows between the leads at a given voltage, 
and at arbitrary asymmetry of the coupling to the two 
leads. We find that the value of the current also allows 
one to distinguish between different stable superconduct- 
ing states at a given voltage. 

The time-dependent calculations revealed that once 
the bias voltage becomes constant in time, the system 
always relaxes into one of the stationary states. Non- 
stationary behavior of physical quantities always decays 
in time, unlike in the case of a DC-biased Josephson junc- 
tion. (Despite it being a non-linear system, a supercon- 
ductor driven by a voltage is therefore fundamentally dif- 
ferent from a nonlinear pendulum driven by an external 
force.) If the bias voltage is changed slowly, an initial 
stationary state evolves adiabatically. By changing the 
voltage slowly we have observed the expected hysteresis 
associated with the existence of two stable states at some 
voltages. 

The rest of the paper is structured as follows. In Sec.HTl 
we specify the model to be studied, and present the equa- 
tions that determine its state. In Sec. IIIII we solve these 
equations analytically, assuming that the system is in a 
stationary state. We analyze the stationary stationary 
states we find and calculate the I-V characteristic of the 
system. In Sec. II VI we establish that the stationary states 
are the only stable states of the DC-biased system. We 
do so by studying the dynamics of the system after a 
perturbation. In Sec. [V] we summarize our main results. 



II. MODEL 

As stated in the introduction, wc consider a supercon- 
ducting island connected to two normal leads by means of 
low transparency tunnel barriers. The superconducting 
order parameter is taken to be spatially isotropic inside 
the island. The physical requirements for this condition 
to hold have already been discussed in the introduction. 
We assume that the dominant energy relaxation mecha- 
nism for the superconductor is tunneling of electrons to 
the leads. For given barrier transparencies, this restricts 
the size of the superconductor to less than the inelastic 
scattering length of quasi-particles inside the supercon- 
ductor. 

Our analysis of the system is based on the Keldysh 
Green function technique i 16 ' 17 ' 18 We start our discussion 
of the equations governing the system by defining the 
necessary Green functions. 
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A. Definition of Green functions 

The Green functions are expectation values of prod- 
ucts of the Heisenberg operators a m ± (t) and a m ± (t) that 
create and annihilate electrons in levels of the isolated 
island. Here m labels single particle levels. The ± in- 
dex accounts for Kramer's degeneracy. As we are dealing 
with a problem involving superconductivity, all Green 



functions are 2x2 matrices in Nambu space. It is useful 
to define Nambu space matrices r/j, j = 0,...,3 such 
that 770 is the identity matrix and 771, r) 2 and rj 3 are 
the standard Pauli matrices. We also define matrices 
V± = (Vi ± W?2)/2. 

The retarded (R), Keldysh (K) and advanced (A) 
Green functions of each level are defined as 1 ^ 



Rm{t, t ) 

K m (t,t') 

Am{t, t ) 



, {a m+ (t),aj n+ (t')} {a m +(t),a m -(t')} 
13 \\{al_(t),al + (t')} {al_(t),a m .(t')} 

(t), a m -(t')] 

l3 \\[al_(t),al + (t')} [al_(t),a m .(t')] 



6(t-t'), 



(2.1a) 

(2.1b) 
(2.1c) 



The Green functions are grouped into a matrix 



B. Equations of motion 



G m (t, t ) 



R m (t,t') K m (t,t') 
A rn {t,t') 



(2.2) 



This further 2x2 matrix structure is referred to as 
Keldysh space. As with Nambu space, it is useful to de- 
fine matrices Tj, j = 0, . . . , 3. The matrix Tj is the same 
as the matrix r\j but now operating in Keldysh space. We 
also carry over the definition of t± from Nambu space. 
A basis for the 4x4 matrices that result from combining 
Keldysh and Nambu indices is constructed by means of 
a tensor product Tj <X> T]k, with the t's always acting in 
Keldysh space and the 77's in Nambu space. 

The quantities that we calculate, namely the order pa- 
rameter A(t) and the current I(t), are collective in the 
sense that they result from the sum of the contributions 
of all the individual levels. Accordingly a formalism ex- 
ists that does not require knowledge of the Green func- 
tions of individual levels but only the sum o 17 ' 20 ' 21 ' 22 ' 23 



l8 s 



Y,Gm{t,t'), Q = G, R, K, A, (2.3) 



that are known as quasi-classical Green functions. Here 
S s is the mean level spacing of the island. 

We will work with the quasi-classical Green functions 
throughout the present section. The advantage of doing 
so is that the theory can be formulated with the least 
amount of clutter. When doing time-dependent numer- 
ics in Sec. IIVI however, we find it more convenient to 
work with the Green functions of the individual levels. 
In principle though, the theory outlined in this section, 
following as it does from the theory outlined in Sec. HV| 
gives exactly the same answers. 



The equations that determine the Green functions can 
be derived from the circuit theory of non-equilibrium 
superconductivit y 21 ' 22 ' 23 Viewed as a matrix in time, 
Nambu and Keldysh indices, the Green function G sat- 
isfies the commutation relation 24 



[H -£,£?] = 0. 



(2.4) 



Here H describes the dynamics of the isolated supercon- 
ductor: 



H(t,t') = t ® m S(t - t') [id t - h(t)} , 



h(t) = 



_ ( -/i,(t) A(t), 



(2.5a) 
(2.5b) 



The matrix h(t) is a remnant of the Bogoliubov-de 
Gennes Hamiltonian. 11 Bearing in mind that we con- 
sider a non-equilibrium setup, we must allow the order 
parameter A(t) and the chemical potential /i s (i) of the 
superconductor to be time-dependent. Their values at 
each instant in time are determined by imposing self- 
consistency. 

The time derivative standing to the right of G in the 
term GH of Eq. (|2.4p can be shifted to act on the second 
time argument of G at the cost of a minus sign, i.e. 

dtG(t,i)d- t 5{t-t') = - j ' did- t G{t,i)8(i-t') 

= -d t >G(t,t'). (2.6) 

The self-energy contains a term corresponding to each 
lead, i.e. 



(2.7) 



I and r referring to the left and right leads respectively. 
The leads act as reservoirs, broadening the island levels 
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to a finite lifetime and determining their filling. The self- 
energy of lead j is = —iTjG^', where Green function 
G<J) of lead j is defined similarly to the Green function 
of the superconductor (Eq. 12.3ft , with the sum now run- 
ning over states in the lead. Here Tj is the tunneling 
rate from any island level to lead j. (For simplicity, we 
take the rates associated with different levels to be the 
same.) The leads are large compared to the supercon- 
ductor, and therefore Gj does not depend on the state 
of the superconductor. Furthermore, since the leads are 
normal, the off-diagonal Nambu space matrix elements 
of the lead Green functions are zero. Explicitly then, the 
Green function for lead j — I, r has the form 

^fcO = <2 ' 8 > 

with 

R U) (t,t') = 5(t-t') m = -A&(t,tf), (2.9a) 

K«\t,t>)=2^f ) aji'tr)- (2 ' 9b) 

The function Oj describes the distribution of particles in 
lead j. In general it is given by 

a j (t,t')= { ^ e -^(*-*')[ 1 _2/,(S)] e -4^W-^(*')] ) 

J 27T 

(2.10) 

where fj (E) is the filling factor of states at energy E in 
lead j. The phase <f>j sets the time-dependent chemical 
potential (ij(t) = d t (j)j(t) in lead j. The time-dependent 
bias voltage between the leads is 

V(t) = (w(f)-M*))/e. (2-11) 

where e is the electron charge. It is convenient to de- 
fine the total inverse lifetime or Thouless energy i?xh = 
Ti + r r and a dimensionless symmetry parameter 7 = 
(Ti — r r )/i?Th- For a perfectly symmetric coupling to 
the leads, 7 = while 7 = ±1 corresponds to the island 
being coupled to only one of the two leads. 

The commutator equation (|2.4[) on its own is not 
enough to specify G uniquely. Indeed what Eq. (|2.4[) 
says is that G has the same eigenstates as H — S, but it 
does not say anything about the eigenvalues of G. Ad- 
ditional to Eq. (|2.4[) there is a also relation between 
the eigenvalues of G and those of H — S. 25 Let |A) be a 
simultaneous eigenstate of H — S and G, such that its 
eigenvalue with respect to H—H is A. Then its eigenvalue 
with respect to G is sgn (Im(A)). (One can show that the 
eigenvalues of H — £ come in complex conjugate pairs 
and that there are no purely real eigenvalues.) Hence G 
squares to unity, i.e. 

G 2 = I. (2.12) 

C. Gauge invariance 

At this point we have defined three different Fermi- 
energies, namely that of the superconductor n s (t) and 



those of the leads [ij (t), j — I, r. Since the reference point 
from which energy is measured is arbitrary, there is some 
redundancy. This redundancy is encoded in a symmetry 
of the equations for the Green function and boils down 
to gauge-invariance. Consider a transformation on the 
Green function 

G^G = UGU^, (2.13a) 
U(t, t') = 5(t - t') t ® exp{i m A(t)). (2.13b) 

As is easy to verify, G obeys equations of the same form 
as G, with chemical potentials and the order parameter 
transformed according to 

H (*) A; it) = H (t) + d t A(t), j = s, I, r, 

(2.14a) 

A(t) -> A(t) = A(t)exp[2iA(t)]. (2.14b) 

When considering stationary solutions we will fix the 
gauge by demanding that A is time-independent. When 
considering non-stationary solutions we will fix the gauge 
such that the reference point from which chemical po- 
tentials are measured is halfway between the chemical 
potentials of the reservoirs, i.e. fJ, r n)(t) = +(— )eV(t)/2. 



D. Self-consistency of A 

The value of the order parameter is set by the self- 
consistency condition 

A(t) = gS s (q m _ {t)a m+ (t)) 

1 n 

= -^Tr^Kfrt)], (2.15) 

where g > is the dimensionless pairing interaction 
strength. This self-consistency equation suffers from the 
usual logarithmic divergence which requires regulariza- 
tion by introducing a large energy cut-off E c . We define 
Ao as the order parameter of an isolated superconductor 
at zero temperature for given g and E c o _. 

Eco 1 f Eca dE 
A = — ^t- =>•-=/ ; 2.16) 

sinhi g J + A 2 , 

This definition then allows us to express A in Eq. (|2.15[) 
in terms of Aq rather than in terms of E c _ , and g. 



E. Current and chemical potential 

The current from the superconductor into reservoir j 
is 2 ^ 



■5 



dt'Tr 



® m (G(t, (t', f) - G w (i, i')G(t', f) 



G K r) = (1 + R7) 



Here Gj is the tunneling conductance of the tunnel bar- 
rier between lead j and the superconductor, and we have 
indicated in square brackets a factor of Ti which equals 
unity in the units we use throughout the paper. The to- 
tal rate of change of the charge in the superconductor 
equals minus the sum of the currents to the leads, i.e. 

-^-Q(t)=Ii(t)+I r (t). (2.18) 
at 

The charge in the superconductor is related to the chem- 
ical potential /i s by means of the capacitance G of the 
superconductor, so that /i s has to obey 



(2.19) 



When the system is not stationary, this equation sets the 
value of /J, s (t) at each instant in time, since dQ(t)/dt can 
be calculated directly from G(t,t'). 



F. Summary 

In summary then, our task is to find the Green function 
G as defined in Eq. (|2.3[) of the superconductor. In gen- 
eral, the procedure for doing this is as follows: We make 
an Ansatz for the order parameter A(i) and the chem- 
ical potential /i s (t). We then diagonalize the operator 
H — £ (that depends on A and fi) . The Green function 
G is constructed in the eigenbasis of H — S, according 
the prescription of Sec. Ill Bl Subsequently we judge the 
correctness of the Ansatz for A(£) and /z(f) by inquiring 
whether Eqs. (|2.15|) and (|2.19j) are satisfied. 



III. STATIONARY SOLUTIONS 

We consider a time-independent bias voltage between 
the left and right reservoirs. In this case the chemical po- 
tentials fii and /i r of the reservoirs are time- independent. 
We make the Ansatz that the chemical potential /i s and 
the order parameter A of the superconductor are also 
time-independent. The Green function G(t,t') only de- 
pends on the time-difference t — t' . It is convenient to 
work with the Fourier transformed Green function G{E) 
which is related to G(t,t') by 

dE 



G(t,t') 



-iE(t-t') 



2tt 



G(E). 



(3.1) 



It is also convenient to construct a traceless operator 
M = H — £ — ro (8 i]o fj, s with Keldysh structure 

M R M k 
M A 



M 



(3.2) 



E T h e 2 

5s [ny 



(2.17) 



I 

In the energy representation the components of M have 
the explicit form 



M R (E) 
M A (E) 
M K (E) 
a(E) 



E + iE Th 
A* 



E - iE Th 
A* 



2iE Th 

1-7 



-A 
E - iE Th 

-A 
-E + iE Th 

a{E) 
a(-E) 

1+7 



-a r (E). 



(3.3a) 
(3.3b) 
(3.3c) 
(3.3d) 



We take the left and right leads to be in local zero- 
temperature equilibrium at Fermi energies /i; = /i + eV/2 
and n r = /i — eV/2 so that the filling factors in both 
reservoirs is a step function fj(E) = 9{—E) and from 
Eq. (|2~TTJ1) follows 



a r (E) 



sgn (E-fi 
sgn (E — fi- 



eV/2) , 
eV/2) , 



(3.4a) 
(3.4b) 



where /i is the average chemical potential (/i r + (ii)/2 
in the leads, in the gauge where the phase of the or- 
der parameter is time- independent. The value of \x will 
later be determined by requiring self-consistency of the 
order parameter A. The Green function G(E) obeys 
[M(E),G(E)\ = 0. The retarded, advanced and Keldysh 
components of this equation are 



[M R (E),R(E)} = [M A (E),A(E)} = 0, 
M R (E)K{E) + M K (E)A(E) 

- R{E)M K {E) - K{E)M A {E) = 0. 



(3.5a) 



(3.5b) 



With the aide of the prescription below Eq. (|2.1ip for 
choosing the eigenvalues of G, one then readily finds for 
the retarded and advanced Green functions 



R{E) 

ME) 
c(E) 



1 



E + iErh 



-A 



c(E) V A* -E- iEr h 

1 / E - iE T h -A 
c(-E) V A* 



(3.6a) 



E + l E T h " (3 ' 6b) 



y/(E + iE Th ) 2 -\A\ 2 . 



(3.6c) 



The function c(E), which we will frequently encounter, 
is defined with branch cuts along the lines E± — ±|A| ± 
x — iETh with x real and positive. The branch with 
\\m.E<=R^±oo C (E)/E = 1 is taken. Considered as a func- 
tion of real E, the real part of c(E) is odd, and the imag- 
inary part is even and positive so that 



2(E)* - -c(-E). 



(3.7) 
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E/\A\ 



FIG. 1: The function c(E) as defined in Eq. (|3.6c|l frequently 
appears in expressions associated with stationary solutions. 
The solid line represents the real part and the dashed line the 
imaginary part. The Thouless energy was taken as _&rh = 
0.1 IAI. 




The real and imaginary parts of c{E) is plotted for real 
E in Fig. □ 

Note that R{E) 2 = A(E) 2 = -q Q as required by 
Eq. 12.121 In general the superconducting density of states 
is g{E) = Trr/ 3 [R(E) - A(E)] /25 s so that we find from 
the solutions for R and A (Eq. I3.G[) 



9(E) 



-Re 



E + iE- 



Th 



c{E) 



(3.8) 



The density of states for an isolated superconductor 
has singularities at energies E — ±|A| of the form 
1 / \/E 2 — |A| 2 . The coupling to the leads regularizes the 
singularities at an energy scale of i?Th- Furthermore, 
whereas the density of states of the isolated supercon- 
ductor vanishes for energies \E\ < |A|, the coupling to 
the leads softens the gap so that there are some states 
for energies \E\ < |A| as shown in Fig. [21 



The next step is to solve Eq. (|3.5bp for K(E). Here 
we use the fact that R{E)K{E)A{E) = -K(E), which 
follows from the requirement that G 2 = I (Eq. I2.12jl . 
Note also that R{E) = M R (E)/c{E) and A(E) = 
Ma(E)/c(—E). Using these identities and multiplying 
Eq. (|375bl) from the left by R(E), we find 



K(E) 



1 



c{E)+c(-E) 
After some algebra we obtain 



[M K {E)-R{E)M K {E)A{E)\. 

(3.9) 



FIG. 2: The density of states of the superconducting island 
(Eq. 13 - 8 P for finite Thouless energy (solid line). The dashed 
line shows the density of states of the isolated superconduct- 
ing island with the same |A|, while the horizontal dot-dashed 
line shows the density of states of the normal island. A value 
of _Etii = 0.1 IAI was used. 



K(E) = 



-KW(E)* 



where 



KW{E) 
K^(-E) 



(3.10) 



K (1 \E) 
K {2) {E) = -ARe 



IAI 2 

S s g{E)a{E) - l -^-Re 



E 



c(E) 



[a(E)+a(-E)} 



c(E) 



[a(E)-a(-E)}- 



iE 



Th 



E 



a(E)+a(-E)} 



(3.11a) 
(3.11b) 



Having obtained K(E) we can find |A| and /x from the 
self-consistency condition Eq. (|2.15p . Below we write the 
real and imaginary parts of the self-consistency equation 



separately. The real part reads 



= dE (Re 



1 



c(E) 



a{E) - a(-E) 



(3.12) 



7 



while the imaginary part reads 



dE^-Re 
E 



1 



c(E) 



[a(E) + a(-E)} 



(3.13) 



These integrals can be done explicitly. We use the iden- 
tities 



dE>Re- ( L- ) =F R (E) 



dE '^yET) 



Fr(0), 
1 



■■Fi(E), 



Th 



where 



Fr(E) = In 



£ + iS T h + c(E) 



Fj(E) = arctan 



A 

Re[c(E)} 



VE; 



Th 



(3.14a) 

(3.14b) 
(3.14c) 

(3.15a) 
(3.15b) 



Here the branch for which —tt/2 < arctan(x) < ir/2 is 
implied. Thus we obtain the transcendental equations 



= (1 - l)F R {^ + M ) + (1 + 7 )F fl ( M - (3.16a) 
. M ) + (l + 7 )F J (/i-^) ) (3.16b) 



= (l- 7 )F>(f 



that determine |A| and /i for given V , £/rh and 7. Below 
we solve these equations analytically in certain limiting 
cases and numerically for more general cases. Only the 
amplitude of A is fixed by these equations. By choosing 
the appropriate gauge [cf. Eq. (|2.14b|) ]. we can set the 
phase of A to any value. In the rest of this section we 
therefore drop the absolute value notation, and take A 
real and positive. 

Before explicitly finding A and /x, we calculate the cur- 
rent from Eq. (|2. 1T|) and the solution for K(E). As- 
suming that the self-consistency equation (Eq. 13 . 13|) is 
fulfilled, we find that the current I r from the supercon- 
ductor to the right lead equals minus the current // from 
the superconductor to the left lead, as it should. For the 
current I = I r = —Ii from the left lead to the right lead 
we find 



/ = 



(l- 7 2 ) eJ E T h 



eV 



Gn 
e 



eJL dE e( E ) 
2 

Re[c(/x+f )-c(/x-f )]. (3.17) 



Here Gjv is the series conductance of the tunneling bar- 
riers to the leads 



Gn — [G ; 1 + G r x ] 



(3.18) 



and Gi and G r are the junction conductances given in 
Eq. f2T7) . 

Now we investigate the transcendental equations (Eqs. 
I3.16[) for /i and A. There are three parameters, namely 



i?Thj 7 and V that determine the solution. Two of these, 
£/rh and 7, are fixed for a given device, while the volt- 
age V can be varied for a given device. (Recall that 
E Th measures the overall coupling to the leads, while 7 
measures the degree of asymmetry between the two lead 
couplings.) Hence it is natural to specify values for E^ 
and 7 and then consider A, /i and / as functions of V. 
In Fig. [3] we show four curves of A versus V, each cor- 
responding to a different choice of the parameters -&rh 
and 7. In Fig. [4] we show the corresponding curves of \x 
versus V. 



1.0 



< 
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1.0 

eV/A 
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2.0 



FIG. 3: The order parameter A versus voltage V, for given 
i?Th and 7. Curves A, B, C and D respectively correspond 
to E Th = 0.35A and 7 = 0.2; E Th = 0.2A and 7 = 0.075; 
E T h = O.IAq and 7 = 0.1; and = 0.01A and 7 = 0.3. 




FIG. 4: The chemical potential fi versus voltage V, for given 
ETh and 7. Curves A, B, C and D correspond to the same 
parameter values as in Fig. [3] 



Let us firstly note the general trend that increasing 
-©rh leads to a smaller order parameter. The reason for 
this is that E^ is the typical time an electron remains 
in the superconductor. The shorter this time (the larger 
.©Th) the harder it is for electrons to form Cooper pairs, 
and superconductivity is inhibited. Secondly, note that 
at large enough i/ph the order parameter is a decreasing 
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function of V. We can therefore obtain the critical Thou- 

(c) 

less energy E T1 J beyond which superconductivity vanishes 
by setting V to zero and asking how large can we make 
£/Th before A becomes zero. 

In the case of V = 0, the self-consistency equations are 
solved by /x = and 

A= { A 0v /l-2^ E Th <A /2 
[ E Th > A /2 

From this we conclude that the critical Thouless energy 

is 4h = Ao/2. 

Having established the range of Exh in which super- 
conductivity persists, we now take a closer look at A as 
a function of V. We have chosen the parameters of the 
four solutions in Fig. [3] to show all the different possible 
shapes that curve of A versus V can take. We see that 
at a given voltage V there can be either zero, one, two 
or three non-zero solutions A. 

To characterize the different types of curve, we con- 
sider V as a function of A on the interval A 6 
[0, A y/1 - 2£ , T h/A ] • In curves of the type A in Fig. El 
V is a monotonically decreasing function of A. In con- 
trast, curves of type B, C and D have local extrema. A 
curve of type B has a local minimum at the left bound- 
ary A = of the A interval on which the function V(A) 
is defined. Then the curve reaches a maximum at some 
intermediate value Ai, before dropping to zero at the 
right boundary A = Aq^I — 2Exh/Ao. Curves C and 
D are distinguished from curve B by the fact that V 
reaches a local maximum instead of a minimum at the 
left boundary A = of the A interval. There is another 
local maximum at intermediate Ai before V drops to 
zero at A = Ao \/l — 2Exh/Ao- In curves of type C, the 
absolute maximum of V as a function of A is at the inter- 
mediate value Ai while for curves of type D the absolute 
maximum of V is at A = 0. 

Next we ask how the Ext, — 7 parameter space is di- 
vided into regions A, B, C and D corresponding to the 
respective types of solution of the self-consistency equa- 
tions. Specifically, which regions share a mutual border? 
Assuming that the function V(A) changes smoothly as 
Exh and 7 are varied, the transitions A <-> B, B «-> C , 
C <-> D and D <-> A are possible. The transition A «-> C 
is not possible. Whenever one tries to smoothly deform 
a curve of type A in Fig. [3] to a curve of type C, one 
invariably reaches a curve of type B or D during an in- 
termediate stage of the deformation. Similarly the tran- 
sition B <-> D is impossible. A smooth deformation of 
a curve of type B into one of type D passes through an 
intermediate stage where the curve is of types A or C. To 
illustrate these ideas we consider a polynomial equation 
of the form 

A„~ Ao 6 Uoy 4 UcJ 2 UcJ ' ( ' 



A, where V(A) is of type A, is given by a > 0, b > 
or a < 0, b > a 2 /4. Region B where V"(A) is of type 
B consists of all points (a, b) such that b < 0. Region 
C consists of all points (a, b) such that a < and < 
b < 3a 2 /16. Region D consists of all points (a, b) such 
that a < and 3a 2 /16 < b < a 2 /A. The regions and 
their borders are shown in the inset in Fig. [5] The most 
pertinent feature of the figure is that the four distinct 
regions meet in the single point a = b = 0. 

Based on a combination of numerical and analytical 
results we have concluded that the Exh — 7 parameter 
space has a very similar topology to this polynomial ex- 
ample. (In principle it could have differed from the poly- 
nomial example by having disconnected regions of the 
same type, for instance two islands of region D, one em- 
bedded in a sea of region A, the other in a sea of region 

C. ) Fig. [5] is a schematic diagram of how the -&rh — 7 
parameter space is partitioned into regions A, B, C and 

D. The following features of the diagram are conjectures 
based on numerical evidence: (1) The regions of types 
A, B, C and D are simply connected. (2) The border 
between regions A and D starts at the corner 7=1, 
■Etii = 0. Other features are deduced from analytical 
results: (1) The line 7 = 0, E Th < A /2 v / 2 belongs to 
region B. (2) The line 7 = 0, A /2V^ < E Th < A /2 
belongs to region A. (3) For E^ > Ao/2 the system 
is in the normal state while it is superconducting for 
-Exh < Aq/2. 4) The border of regions D and C meets 
the border of region B and C at Exh = 0, 7 = 0. 
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FIG. 5: Schematic diagram of the partitioning of the Exh — 7 
parameter space into regions where the curve of A versus V 
is of the types A, B, C and D (Fig.[3]|. The regions A, B, C, 
and D meet at point Q. The line Exh = Ao/2 separates the 
normal and superconducting regions of parameter space. The 
dots in the figure indicate the parameter values that corre- 
spond to the curves in Fig. [3] The inset shows the regions A, 
B, C, and D in the parameter space a — b of the polynomial 
V7A) of Eq. (pT2"Uj) . The topology of the E T h — 7 parameter 
space of the superconductor in the region of the point Q can 
be understood by considering the topology of the parameter 
space of the polynomial. 



We ask what are the respective regions of the a — b plane 
in which V(A) is a curve of type A, 5, C and D. Region 



In region D, superconductivity can persist up to volt- 
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ages that are large compared to Ao. For given i?xh and 
7 there is however always a critical voltage V c beyond 
which superconductivity ceases. (This is a second order 
phase-transition.) For the voltage V c we have obtained 
the following analytical result from Eq. (|3 . 16[) . At finite 
7 and for 2?Th sufficiently small, V c obeys the power law 



y{c) 



Ao 
2e 



2E Th ttA 
— — sec — — 
A 2 



A = 



l + M" 



(3.21) 



This power law is valid as long as ^> Ao/e. It is 
from this result that we are able to conclude that the 
region of finite 7 and infinitesimal i?xh belongs to region 
D. 

Another analytical result can be obtained for the case 
of perfectly symmetric coupling to the leads, i.e. 7 = 0. 
In this case, Eq. (|3.16bj) is solved by n = and the 
relation between A and V can be stated as 



eV = Ao ( 1 + 



0/ \ 



1 - 



A 2 



(3.22) 



This result is plotted for several values of Et\i in Fig. [6] 
It is from this result that we are able to conclude that the 
line segment 7 = 0, < i?Th < Ao/2\/2 belongs to region 
B while the line segment 7 = 0, A /2v / 2 < E T h < A /2 
belongs to region A. The E Th — ► limit of Eq. (|3.22|) can 
be obtained by considering a bulk superconductor and 
assuming a quasi-particle distribution function n(E) = 
(6{-eV/2-E)+6{eV/2-E))/2. It is also worth noting 
that the same result is obtained for a T-junction where 
the stem of the T is a superconductor and the bar a 
voltage-biased dirty normal metal wire.— 

Finally, we consider the I — V curves associated with 
the solutions A and fi of Figs. [3]and|3J The results are 
shown in Fig. [7] From these curves we can infer the re- 
sults that will be obtained in an experiment in which 
the voltage V is swept adiabatically from zero to several 
Ao/e and back to zero. In region A of parameter space 
there is a single current associated with each voltage. At 
some voltage V+ of order Ao/e the system makes a phase 
transition to the normal state, but this does not lead to 
a discontinuity in the current versus voltage curve. In 
contrast, in regions B, C and D, the current will make 
discontinuous jumps as the voltage is swept. Hysteresis 
will also be observed. Suppose the device is in region B 
of parameter space. As V is swept from upwards, a 
voltage V+ is crossed where the current makes a finite 
jump. After the jump, the system is in the normal state 
and the current is GjyV . (Gn is the normal state con- 
ductance of the setup, cf. Eq. 13.181 ) On the backward 
sweep from V > V+ to zero, the system remains normal 
when V+ is reached. At some voltage V- < V+ the cur- 
rent jumps from its value GnV— in the normal state to 
a smaller value, signaling the onset of superconductiv- 
ity. The behavior of the system in region C of parameter 
space is similar. The upward sweep of the voltage pro- 
duces a jump in the current at a voltage V+. After the 




FIG. 6: The order parameter A versus voltage V for symmet- 
ric coupling to the leads, i.e. 7 = 0, according to Eq. (|3.22|l . 
Different curves correspond to different Exh- From the top 
curve to bottom curve we took i?Th/Ao = .01, 0.14, 0.26, 
l/2 v / 2(^ 0.35), 0.42 and 0.47. The curve corresponding to 
I/2V2A0 is plotted thicker than the others. For 
smaller _&rh are of type B with two non-zero values for A 
at some voltages. For larger i?Th, curves are of type A, with 
at most one non-zero A at every voltage. 



jump the system is normal and the current is given by 
/ = GnV. The difference from region B appears when 
the voltage is swept back from V+ to zero. At some volt- 
age smaller than V+ the current starts deviating from its 
value in the normal state, but there is no discontinuous 
jump yet. Even so, the system has turned superconduct- 
ing. When the jump in current now occurs at V- < V+, 
the system switches between two different superconduct- 
ing states. Finally, for parameters in region D, the volt- 
age sweep produces results similar to that in region C. 
The difference between regions C and D is that in D the 
system also jumps between two superconducting states 
at V+ during the forward sweep. 



IV. DYNAMICS 

We concluded the previous section with a discussion 
of hysteresis in the current -voltage characteristic of the 
superconducting island. The conclusions we drew rely on 
the assumption that after the system is perturbed by a 
change in the bias voltage, it relaxes into a stationary 
state. The validity of this assumption is by no means ob- 
vious, since the system is driven (by the bias voltage) and 
the stationary state is not an equilibrium state. Frankly, 
our own initial expectation was that the presence of a 
bias voltage would cause the dynamics of |A(t)| to be 
quasi-periodic or chaotic. We therefore did numerical 
simulations in order to investigate the dynamics of |A(f)| 
in the presence of a bias voltage. Our main result is this: 
Suppose the bias voltage assumes the constant value V/ 
for times t > tf. Then (contrary to our original expecta- 
tions) at t ^> tf the superconductor will always be found 
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FIG. 7: The current / through the superconductor versus 
the voltage V across it. Curves A, B, C and D correspond 
to the respective parameter values quoted in Figs. [3] and [4] 
The dashed line shows the current through the system in the 
absence of superconductivity. 



in one of the stationary states associated with Vf . This 
is true regardless of the history of the system prior to 
t < tf. In particular, the time dependence of the bias 
voltage prior to tt does not matter. Nor does the state 
of the superconductor prior to tf matter. Only when 
there is more than one non-zero stationary solution asso- 
ciated with Vf does the history of the system have any 
baring on its final state. In this case, the history of the 
system determines which of the possible stationary states 
eventually becomes the final state of the superconductor. 
For slowly varying voltages, the predictions of the previ- 
ous section regarding hysteresis are confirmed. In this 
section we discuss the numerics that yielded the above 
results. 

For the purpose of numerics we find it advantageous 
not to take the sum over levels of the Green function as 
we did in the previous sections. Instead we work with the 
Green functions of each individual level. The advantage 
of this scheme is that it allows us to work with ordinary 
differential equations. From these differential equations 
it is straight-forward to construct a time-series in which 
the next element can be calculated if the present elements 
are known. As far as we can see, no such 'local in time' 
update equations exist for the Green functions summed 
over levels. Naturally there are disadvantages to work- 
ing with the individual level Green functions as well; the 
number of equations to be solved numerically is increased 
significantly. As a result the calculation is computation- 
ally expensive and therefore time-consuming. 

The Green functions of the individual levels obey the 
equations 

(H m - E)G m = G m {H m - S) = I. (4.1) 

Here H m differs from the operator H that appeared in 
Eq. (|2.5[) in that it contains the energy e m of level m. It 



is explicitly given by 

H m (t,t') = T ® m S{t-t')[id t -h m (t)}, (4.2a) 

The operator h m (t) is the time-dependent Bogoliubov-de 
Gennes Hamiltonianfii The self-energy E is the same as 
in Sect. EE 

We measure energies from a point halfway between 
the chemical potentials of the leads. As a result the 
phases (t) that appear in the reservoir self-energies are 
0r(o(*) = ~K — )<?K*)/2 where (f> is related to the voltage 
V by V(t) = d t <P{t)/e. 

We parameterize the Green functions in terms of a set 
of auxiliary functions. This eliminates some redundan- 
cies that are present due to symmetries of the equations 
of motion. We start by noting that since the retarded 
and advanced Green functions are related by Eq. (|2.1cD . 
we do not need to consider both. We work with the re- 
tarded Green function. We define a matrix r m (£, r) that 
is related to R m {t, t — r) by the equation 

R m (t,t-T) = ir) 3 r m (t,r), (4.3a) 
r m (t,r) = r$(t,T)rio+ir m (t,T).Ti. (4.3b) 

Here the component (<, r) of r m (t, r) is a scalar func- 
tion whereas the other three components are grouped into 
a vector r m (t, r) such that 

r m (t,r) = (rW(t,r),rW(t,r),rW(t,T)) . (4.4) 

The vector r] — (771 , 772 , 773 ) contains the Pauli matrices in 
Nambu space. Before the voltage bias between the leads 
is established, (i.e. for t < and all r), the functions 

rm{t,r) and r m (t, r) are real. When the equations of 
motion (|4.1[) for the retarded Green function are rewrit- 
ten in terms and r, we find that their reality is pre- 
served at all times. 

Next we consider the Keldysh Green function. In order 
to calculate the time-evolution of the order parameter 
we only need to know the Keldysh Green function at 
coinciding times. Here the parameterization 

K m (t,t) = ir) 3 k m {t)-rj, (4.5) 

in terms of a real vector 

k m (t)=(k£\t),k£\t),kg>(tj), (4.6) 

is respected by the initial condition and preserved by the 
equations of motion. 

From the equations of motion (Eq. 14. ip we derive dif- 
ferential equations 

-jrr m (t,T) = b m (t)r m (t,r) - r m (t,T)b m (t - r), (4.7) 
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—k m {t) + 26 m (i) x k m (t) + 2E Th k m {t) = 4E Th f m (t), 
at 

(4.8) 

for the matrix r m (t,T) and the vector k m (t). The equa- 
tion (O with E Th = was studied in Refs. HH. In 
these references the dynamics of the order parameter of 
an isolated superconductor was calculated. We see that 
coupling the system to leads introduces two terms pro- 
portional to i?Th- One (on the left-hand side of Eq. (14. 8|) ) 
can be considered a damping term and is proportional to 
k m (t). The other (on the right-hand side of Eq. I|4.8p) 



can be considered a driving or source term. 

In Eq. (|4.8p . b m (t) is a matrix and b m (t) a vector such 
that 

b m (t) = ib m (t) ■ r], (4.9a) 
b m {t) = (RcA(i),-ImA(t),/i s (*)-e m ). (4.9b) 

The equation for k m (t) contains a source term 
4£ , Th/ m (i)- The vector / m (t) is given by 



fm(t) 



dr 



»(t, 



,(o) 



(*,t) 



!(t,r) x s(t,r) 



(4.10) 



In this equation the scalar function s' ^ (t, r) and the vec- 
tor s(t, t) parameterize the Kcldysh component of the 
self-energy as follows 

Sif(M') = 2E Th m s(t,r), (4.11a) 

s(t,r) = S {0) {t,T)T] +is(t 7 T) J]. (4.11b) 

Referring back to Sec. [HI where Yik is expressed in terms 
of the Fourier transform of the reservoir filling factors, we 
find explicitly 



a C°>(i,T) 



1 fl\ 6(t) - 6(t - t) , 

- V I- cos *LL-|! 1, (4,12a) 



s(f , T) _ -Jlfro,*. *"-?' ~ T) 

7TT V 2 



(4.12b) 



By imposing self-consistency, the order parameter A(i) 
is expressed in terms of the components of k m (t) as 



(4.13) 



m— — n 



where the number of levels on the island is 2S1 + 1. This 
makes the differential equations non-linear, since they 
contain terms in which A(t) multiplies k m and r. We 
eliminate the dimensionless pairing strength g and the 
mean level spacing S s in favor of A eq , the order parameter 
of the island in equilibrium, by means of the equilibrium 
self-consistency relation 



2 A 1 2 £ 
— — = y arctan 

9 s ' m t^n ^ n 



Em 



where 



A 2 . 



(4.14) 



(4.15) 



As with A(t), the chemical potential ^ s (t) is deter- 
mined by a self-consistency equation. The chemical po- 
tential takes into account the work that must be per- 
formed against the electric field of the excess charge 
on the superconductor in order to add more charge. 
Thus Hs(t) is related to the charge of the island by 
fJ, s (t) — e[Q(t) — Qo]/C where C is the capacitance of 
the island. In this equation Qq represents the fixed posi- 
tive background charge and Q(t) is the combined charge 
of all the electrons on the island. Since the differential 
equations (|4.7p and (|4.8[) only depend on the difference 
Us (t) — /i s (t — r) , the positive background charge need not 
be specified. The charge Q(t) is related to the Keldysh 
Green function. Indeed, the average number n m (t) of 
electrons (with spin-degeneracy included) in level m at 
time t is given by n m (t) = (1 — iTr[K m (t,t)])/2, Hence 
l^ s (t) is related to k m by the equation 



n 



^(t)-^(t-r) = — J2 k£>(t)-k£>(t-r). (4.16) 

rn— — f2 

Finally, we have to specify the initial conditions for 
r m (t,r) and k m (t). We will assume for our simulations 
that the voltage between the reservoirs is zero and the 
system is in zero-temperature equilibrium for times t < 0. 
The corresponding initial condition at t = is 



^ 0) (0,r) = -0(r)e- BT ^cos(e m r), (4.17a) 

rn ) , 

(4.17b) 



r m (0,r) = 9(r)e- E ^^f^(-A cq ,0,s m ), 



MQ) 



arctan m (A cq , 0, — i 

4m 7T Erh 



■m. 



. (4.17c) 



We are now ready to study the time-evolution of A(t) 
when a non-zero bias voltage V(t) between the leads is 
present for times t > 0. In the calculations we report on 
here, we worked with i?Th — 0.069 A and three different 
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FIG. 8: Top panel: The stationary solutions for A vs. the bias 
voltage V, corresponding to the parameters used in generating 
Fig. The outer curve (red) corresponds to E?h = 0.069Ao 
and 7 = 0.2. The middle curve (black) corresponds to Et\i = 
0.069Ao and 7 = 0.1. The inner curve (blue) corresponds to 
Eti, — 0.069Ao and 7 = 0.05. The vertical lines indicate Vf 
and V2. (Vi is beyond the left edge of the figure.) Bottom 
panel: The time-dependence of the voltage. The upper (blue) 
curve corresponds to the blue curves of A vs. t in Fig. [9] The 
lower (red) curve corresponds to the red curves of A vs. t in 
Fig. [9] 



7, namely 7 = 0.05, 7 = 0.1 and 7 — 0.2. These all 
correspond to points from regions C and D in the £/rh - 7 
parameter space of Fig. [5l Hence, for each of the pa- 
rameter choices, there is a bias voltage interval [V_,V+] 
where there are more than one non-zero stationary solu- 
tions for |A|. The three curves of stationary |A| versus 
V, corresponding to the different parameter choices, are 
plotted in the top panel of Fig. [5J 

For given £/rh and 7 we did two numerical runs with 
different time-dependent voltages V{t). The two voltages 
are plotted as functions of time in the bottom panel of 
Fig. [5J In the first run we start by rapidly establishing 
a bias voltage V\ < V-. Rapid here means dV/dt 3> 
Ao-E/rh/e- In this case V changes by an amount of order 
Ao/e — the scale at which the stationary solution for |A| 
depends on V — in a time that is short compared to the 
relaxation time E^£. (Slow refers to the opposite limit.) 
We then keep the voltage constant at V\ for a length of 
time of several -Erh- This time- interval is long enough 
for any transient behavior induced by the rapid change 
of V(t) to disappear. We then slowly increase the bias 
voltage until we reach a bias voltage Vf G for 
which more than one non-zero stationary solutions exist. 
In the second run we start by rapidly establishing a bias 
voltage V2 > V+. We keep the voltage fixed at V2 for a 



time of several E Th . We then slowly decrease the voltage 
to Vf. The values of Vf, V2 and Vf were chosen V\ — 
0.83A , V 2 = 1.76A and Vf = 1.34A . The calculations 
were performed with 501 equally spaced levels with level 
spacing S s = 0.018A and the capacitance was chosen 



C* = 0.1e7A . 
1.0 



0.5 



0.0 



15 



30 



45 




tB 



Th 



FIG. 9: The amplitude of the order parameter as a function of 
time. All curves are for _Exh = 0.069Ao. The top, middle and 
bottom panels correspond to 7 = 0.05, 7 = 0.1 and 7 = 0.2 
respectively. The red curves correspond to a voltage that is 
increased from Vi = 0.83Ao to Vf = 1.34Ao. The blue curves 
correspond to the voltage being decreased from V% = 1.76Ao 
to Vf — 1.34Ao. The vertical lines indicate the time-interval 
in which the voltage changes from either Vi or Vi to Vf . The 
thin horizontal lines correspond to the stationary values of 
A I for a bias voltage V = Vf as calculated from Eq. (|3.16[1 . 



The resulting |A| are plotted as functions of time in 
Fig. [5J They firstly show that after the initial rapid 
change in the bias voltage the system always relaxes into 
a stationary state consistent with the new voltage. The 
relaxation takes a time of the order E^. Secondly, if 
the system is in a stationary state, and the bias voltage 
is changed slowly then | A(i)| adiabatically tracks the sta- 
tionary solution corresponding to the instantaneous value 
of the voltage. This is seen most clearly in Fig. [TU] where 
we plot |A(i)| as a function of V(t) and compare this to 
the stationary |A| vs. constant V curves. Our prediction 
about hysteresis is confirmed. Systems with different his- 
tories end up in different stationary states at the same 
voltage bias. If the voltage is slowly swept from a small 
initial voltage to Vf £ [V-, V+] a stationary state with a 
large value for | A| is reached. If the voltage is swept from 
a large initial voltage to Vf € [V-, V + ], a stationary state 
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is reached that corresponds to a small value of |A|. We 
must mention here that we observe some slow drift (too 
slow to be visible in Fig. [9]) in |A| after the voltage has 
reached Vf. The value of |A| seems to increase linearly 
at a rate d\A\/dt ~ 10 _4 Aq. Within the numerical accu- 
racy of the calculation, this is negligible and we believe 
the drift is simply an artifact of the numerics. 

In our data there is one exception to the rule of adi- 
abatic evolution. In the middle panel of Fig. [9j |A(t)| 
takes much longer that EJ™. to respond when the volt- 
age is changed from Vi to Vf. Hence, |A(f)| as a func- 
tion of V(t) does not track the stationary solution in 
this instance. The reason is the following: For a voltage 
V = V2, the only stationary solution has A = 0. When 
the voltage is decreased to V/, a non-zero stationary so- 
lution for A exists. However A = is still a valid state, 
albeit unstable. The time it takes the system to diverges 
from the unstable state is not determined by E-rh but 
rather by small numerical errors that perturb the unsta- 
ble state. 
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FIG. 10: The amplitude of the order parameter |A(t)| as a 
function of voltage V(t). The parameter values of the three 
panels are the same as those in Fig. [9] i.e. all curves are 
for iJxh = 0.069Ao. The top, middle and bottom panels 
correspond to 7 = 0.05, 7 = 0.1 and 7 = 0.2 respectively. 
The red curves correspond to a voltage that is increased from 
V\ = 0.83Ao to Vf = 1.34Ao. The blue curves correspond 
to the voltage being decreased from V2 = 1.76Ao to Vf = 
1.34Aq. The dashed lines represents the stationary value of 
A vs. V, as calculated in Sec. Mil and plotted in Fig. [8] 



One possible explanation for the observed stability of 
the stationary states is overdamping. According to this 
hypothesis, if we decrease the Thouless energy further, 
thereby decreasing the damping, the stationary solutions 
will become unstable. Some evidence for the hypothesis 
might be visible in Fig. [5] After the voltage is changed 
rapidly, we might expect | A(i)| to perform damped oscil- 
lations while relaxing to the new stationary state. How- 
ever in Fig. [9] no such oscillations are visible, apparently 
implying that the relaxation rate is larger than the os- 
cillation frequency. There is however another possible 
explanation for the lack of oscillatory behavior after an 
abrupt change in V . The argument is that an abrupt 
change in V cannot be communicated to the system 
abruptly, but only at a rate comparable to the damp- 
ing rate Erh- This is because the superconductor learns 
of the change in voltage by the same mechanism as by 
which damping occurs, that is, by tunneling of particles 
between the leads and the island. Hence the response of 
the order parameter is always gradual. 

How do we test whether overdamping hypothesis is 
true or false? Ideally we would have liked to repeat the 
above numerical calculation with a smaller value of Exh 
and see if the stationary states are still stable. However, 
the value Exh = 0.069Ao that we used above is close to 
the smallest value for which we can do reliable numerics 
in reasonable time. Since we cannot make Exh smaller, 
we resolve the issue of overdamping as follows. We com- 
pare the dynamics of A after an abrupt change in the 
pairing interaction strength g at Exh = 0.069Ao to the 
dynamics after a change in g at Exh = 0j2£ We know 
that in the isolated system, (Exh = 0) |A| will perform 
persistent oscillations^ The period of oscillation gives 
a typical time-scale for the internal dynamics of A. If, 
in the open system (i.e. Exh 7^ 0), we observe a few 
damped oscillations (the more the better) in |A| before 
the system relaxes to equilibrium, it means that damping 
occurs at a timescale larger than that of the internal dy- 
namics of the superconductor. In this case the hypothesis 
of overdamping is discredited. 

In our numerical implementation of the above, we work 
with the following parameters: The initial pairing inter- 
action is such that for t < 0, A = A*. The increased 
pairing interaction strength corresponds to an equilib- 
rium value of the order parameter A/ = 20 A;. The per- 
sistent oscillations of |A(t)| in the isolated system are 
shown in the blue curve in Fig. [TTJ We repeat the cal- 
culation, now for a superconductor connected to leads. 
We use a s energy E Th = 0.075A/. The result for |A(t)| 
in the presence of leads is the red curve in Fig. [TTJ We 
see that |A(i)| eventually decays to a constant, as ex- 
pected. The extent of the damping is such that several 
oscillations are completed within the decay time. Hence 
we conclude that the numerical results that we obtained 
previously are outside the regime of overdamping. It fol- 
lows that the lack of oscillatory behavior in Fig. [5] is due 
to the fact that the superconductor only gradually be- 
comes aware of a change in the voltage. 
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FIG. 11: The order parameter versus time after the pairing 
strength was increased from Aj = 0.05A/ to Af abruptly 
at t = 0. The blue curve is for an isolated superconductor 
while the red curve is for a superconductor connected to leads. 
For this case a Thouless energy £xh = 0.075A/ was used. 
The data was obtained using 501 equally spaced levels with 
level spacing 5 S = 0.02 A/. The capacitance was chosen C = 
0.1e a /A/. 



V. CONCLUSION 

We have studied a voltage biased NISIN junction i.e. 
a superconducting island connected to normal leads by 
means of tunnel junctions. We restricted ourselves to the 
regime where the dominant energy relaxation mechanism 
in the superconductor is the tunneling of electrons from 
the superconductor to the leads. We also restricted our- 
selves to the regime of low transparency junctions where 
the position dependence of the order parameter inside 



the superconductor can be neglected. 

In Sec. lIIII we found the stationary states of the system. 
For these, the order parameter A and the chemical po- 
tential are implicitly determined by Eq. (|3.16p . We also 
found the current between the leads [cf. Eq. (|3.17p ] . The 
most striking feature of the stationary states is that there 
can be more than one stationary state at a given voltage. 
These are characterized by different values of |A| and of 
the current as can be seen in the I-V curves of Fig. [7] 
Depending on system parameters, superconductivity can 
survive up to voltages large compared to A , the order 
parameter of the isolated superconductor. In this case, 
increasing the voltage eventually leads to a second order 
phase transition to the normal state. We have found that 
the critical voltage at which the transition occurs obeys 
a power-law [cf. Eq. (|3.21[) ] . 

In Sec. [TV] we studied time-dependent states of the 
system. In this way we were able to demonstrate the 
stability of the stationary states we have found in the 
previous section. Our results also indicate that a DC bi- 
ased system always relaxes into a stationary state. In 
the parameter region of multiple stationary states we 
demonstrated bi-stability. Associated with this are first 
order phase-transitions: there are critical voltages where 
A (and the current) make finite jumps. Furthermore, 
there is hysteresis of |A| and the current associated with 
the bi-stability. 
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